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ABSTRACT 

We present numerical computations of the equilibrium configurations of tidally-locked homogeneous 
binaries, rotating in circular orbits. Unlike the classical Roche approximations, we self-consistently 
account for the tidal and rotational deformations of both components, and relax the assumptions 
of ellipsoidal configurations and Keplerian rotation. We find numerical solutions for mass ratios q 
between 10 -3 and 1, starting at a small angular velocity for which tidal and rotational deformations 
are small, and following a sequence of increasing angular velocities. Each series terminates at an 
appropriate "Roche limit" , above which no equilibrium solution can be found. Even though the 
Roche limit is crossed before the "Roche lobe" is filled, any further increase in the angular velocity 
will result in mass-loss. For close, comparable-mass binaries, we find that local deviations from 
ellipsoidal forms may be as large as 10 — 20%, and departures from Keplerian rotation are significant. 
We compute the light curves that arise from our equilibrium configurations, assuming their distance is 
^ 1 AU (e.g. in the Kuiper Belt). We consider both backscatter (proportional to the projected area) 
and diffuse (Lambert) reflections. Backscatter reflection always yields two minima of equal depths. 
Diffuse reflection, which is sensitive to the surface curvature, generally gives rise to unequal minima. 
We find detectable intensity differences of up to 10% between our light curves and those arising from 
the Roche approximations. Finally, we apply our models to Kuiper Belt binary 2001 QG298, and find 
a nearly edge-on binary with a mass ratio q — 0.931q oL angular velocity cj 2 /Gp = 0.333 ± 0.001 
(statistical errors only), and pure diffuse reflection. For the observed period of 2001 QG298, these 
parameters imply a bulk density, p = 0.72 ± 0.04 g cm -3 . 

Subject headings: Kuiper Belt - minor planets, asteroids - solar system: general 



1. INTRODUCTION 

The Kuiper Belt consists of a large number of small 
objects in heliocentric orbits beyond Neptune. The ex- 
istence of the Kuiper Belt was suggested on theoretical 
grounds by Edgeworth (1943) and Kuiper (1951), but 
it was not until 1992 that the first Kuiper Belt object 
(KBO) was detected (Jewitt & Luu 1993). To date, more 
than 1000 KBOs are known. They are thought to be 
relics of the Sun's accretion disk, and to hold signatures 
of the planetary migration process. Their physical prop- 
erties, including their mass-distribution, compositions, 
and binary-fraction may thus hold the key to our under- 
standing of the formation and evolution of the early solar 
system. A comprehensive review is provided by Luu & 
Jewitt (2002). 

The properties of Kuiper belt binaries place impor- 
tant constraints on theories of solar system evolution. 
In particular, the distributions of separations and mass 
ratios are a unique signature of the binary formation 
process. Over the past decade, more than 50 resolved 
Kuiper belt binaries have been discovered and stud- 
ied (Noll et al. 2007, 2008). These can be broadly di- 
vided into two groups (Noll et al. 2008). The first con- 
sists of small satellites in orbits about larger KBOs. 
These systems are believed to have formed through a 
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large collision, followed by tidal evolution, much like the 
Earth-Moon and Pluto-Charon systems (Hartmann & 
Davis 1975; Cameron & Ward 1976; McKinnon 1989). 
The second - more abundant - group, consists of com- 
parable mass, large separation binaries. These are in- 
consistent with the classical formation mechanisms, and 
have prompted an abundance of new theoretical models 
(Weidenschilling 2002; Goldreich et al. 2002; Funato et 
al. 2004; Astakov et al. 2005; Lee et al. 2007; Schlicht- 
ing & Sari 2008; Gamboa Suarez et al. 2010; Naoz et 
al. 2010). Recently, Sheppard & Jewitt (2004) postu- 
lated the existence of a third group of Kuiper belt bina- 
ries. By following the light curves of KBOs, they were 
able to identify an unresolved, comparable-mass small- 
separation binary-candidate, 2001 QB298, based on its 
variability, photometric range, and period. Their anal- 
ysis indicates that at least 10 — 20% of all large KBOs 
may in fact be unresolved close-binary systems. 

Our current interpretation of the observed light curves 
of KBOs relies on the classical theory of the equilibrium 
figures of rotation. The equilibrium configurations of ro- 
tating fluid bodies is a classical problem that was first 
investigated by Newton in the context of the figure of 
the Earth. Newton showed that the slow rotation of the 
Earth causes it to become a slightly oblate spheroid. In 
1742, Maclaurin generalized the problem, by relaxing the 
assumption of slow rotation. He derived a general rela- 
tion between the velocity of rotation and the eccentricity 
of the resulting equilibrium spheroid. Maclaurin's re- 
lation implies that for any given angular velocity, two 
equilibrium solutions may be found, with different ec- 
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ccntricities. Decades later, Jacobi (1834) realized that 
triaxial ellipsoids are also possible equilibrium configura- 
tions of homogeneous rotating masses. 

The study of possible equilibrium configurations of bi- 
nary systems was introduced by Roche (1847). Roche 
considered the equilibrium form of a satellite rotating 
about a rigid sphere in a circular Keplerian orbit. The 
satellite is deformed both by rotation and by the tides 
exerted by its spherical companion. A more accurate ap- 
proach was later taken by Darwin (1906), who revisited 
the problem, allowing for mutual deformations of both 
bodies. Solutions to Darwin's problem may be found 
only in the limit of equal mass ratio (q = 1) or extreme 
mass ratio (q <C 1). Both the Roche and Darwin solu- 
tions rely on the assumption that the equilibrium con- 
figurations are triaxial ellipsoids. A thorough review, 
analytic derivation, and stability analysis of these clas- 
sical equilibrium configurations can be found in Chan- 
drasekhar (1969). 

In modern astrophysical research, these classical Roche 
ellipsoidal approximations are used in the study of a wide 
range of contact-binaries, including binaries in the so- 
lar system. Weidenschilling (1980) used equal-mass Dar- 
win models to fit the observed light curves of asteroids 
624 Hektor and 216 Kleopatra. In 1984, Leone et al. an- 
alytically constructed arbitrary mass-ratio Roche ellip- 
soidal models. When computing a Roche binary model 
for an arbitrary mass ratio, each component is in turn as- 
sumed to be spherical while the ellipsoidal configuration 
of its companion is calculated. Leone et al. (1984) used 
their theoretical photometric-range as a function of the 
angular velocity to constrain the mass ratios and densi- 
ties of observed variable asteroids. Ccllino et al. (1985) 
explicitly computed the light curves arising from the 
equilibrium Roche models of Leone et al., and compared 
them with an observed sample of asteroid light curves to 
constrain their orbital parameters and densities. Jewitt 
& Sheppard (2002) used the observed period and pho- 
tometric range of KBO 20000 Varuna to consider pos- 
sible (single) Jacobi ellipsoid and Roche binary mod- 
els for this object. Takahashi & Ip (2004) computed 
the light curves arising from Roche ellipsoidal configura- 
tions, to confirm the nature of the suspected binary KBO 
2001 QG298- Lacerda & Jewitt (2007) constructed a li- 
brary of Jacobi-ellipsoid and Roche-binary light curves, 
and investigated the nature and densities of four vari- 
able KBOs (20000 Varuna, 2003 EL 6 i, 2001 QG 298 , and 
2000 GN m ) and of the Trojan asteroid 624 Hektor. 
Descamps (2008) applied Roche binary models to a sam- 
ple of variable asteroids. 

All of the binary models mentioned above rely on the 
assumptions of the Roche approximation, namely that 
each component is deformed by the tides of a spheri- 
cal companion, and that the resulting configurations are 
ellipsoids rotating in Keplerian orbits. In this paper 
we present new numerical computations of the equilib- 
rium configurations of tidally locked homogeneous bi- 
naries, orbiting in circular orbits. We self-consistently 
take into account the tidal and rotational deformations of 
both components, and relax the assumptions of ellipsoidal 
configurations and Keplerian rotation. For comparable- 
mass small-separation binaries, departures from ellip- 
soidal configurations and Keplerian rotation become sig- 
nificant. 



Numerical solutions of the equilibrium configurations 
were previously computed for the cases of general mass 
ratio homogeneous binaries (Hachisu & Eriguchi 1984a) 
and for equal-mass polytropic binaries (Hachisu & 
Eriguchi 1984b). These previous works expanded the 
local gravitational potentials in terms of Legendre poly- 
nomials, which were used to approximate the potential 
on the surface of the two components. In our models 
we explicitly compute the local gravitational and rota- 
tional potentials on the surface of the two bodies, and 
use a Newton-Raphson based scheme to converge to an 
equilibrium solution, for which the surface of each of the 
components is an equipotential surface. We then use 
our numerical solutions to calculate the light-curves that 
these equilibrium configurations exhibit if placed in the 
Kuiper Belt, and demonstrate how these models can be 
used by fitting the observed light curve of Kuiper belt 
binary 2001 QG 298 . 

In our computations we make several simplifying as- 
sumptions. First, we assume that the bodies' gravity 
dominates over their internal strength, so that they take 
the forms of rotating fluid binaries appropriate for their 
angular velocity. The ratio of the material rigidity to 
self-gravity determines a size scale above which bodies 
may be considered gravity-dominated. For rigid mono- 
liths, this size scale is ~ 10 4 km. However, many of the 
smaller solar-system bodies are expected to be "rubble- 
piles" , for which the effective rigidity is reduced (Goldre- 
ich & Sari 2009; and references therein). Rubble piles 
larger than a few hundred kilometers are likely to be 
gravity dominated. In addition, over the age of the solar 
system, even smaller bodies will gradually take the forms 
dictated by gravity. 

Second, we assume that the densities of the two binary 
components are equal. While this is likely the case for bi- 
naries formed through collisions, it is not necessarily true 
for binaries created via three-body interactions or due to 
dynamical friction. At the opposite extreme, when the 
density ratio is infinite, one of the bodies is effectively 
compressed into a point mass while the other takes a 
form similar to the classical Roche solution. These two 
limiting cases therefore span the range of possible con- 
figurations for different density-ratios. 

Third, we assume that the binaries are tidally locked, 
so that their orbital and spin frequencies are equal, and 
their orbits circularized. If the orbital and spin velocities 
differ, the tidal bulge is carried ahead (or lags behind) 
the companion. Torques then act to re-align the sys- 
tem, inducing oscillations about the equilibrium aligned 
configuration. Internal dissipation damps these oscilla- 
tions, and results in synchronized circular orbits. For the 
fiducial physical properties of large KBOs (R <~ 100 km, 
p ~ 2 g cm ~ 3 , fcrubbic ~ 10~ 3 , Q ~ 100), such syn- 
chronization is expected to occur on time scales that 
are shorter than the age of the solar system for bina- 
ries forming with initial separations < 50i? (Goldreich & 
Sari 2009). 

The outline of this paper is as follows. In section 2, 
we briefly describe the numerical method that we use to 
compute the equilibrium configurations as a function of 
the mass ratio and angular velocity. Additional details 
are provided in APPENDIX A. In Section 3 we discuss 
the equilibrium figures of rotation. We first use ana- 
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lytic approximations to clarify the nature of the equilib- 
rium solutions and the evolution surrounding the "Roche 
limit". We then present exact solutions for the equilib- 
rium configurations using our numerical approach. In 
Section 4 we show that our equilibrium configurations 
are generally non-ellipsoidal, and deviate from pure Ke- 
plerian rotation. We explain how the angular momentum 
depends on the mass ratio, angular velocity and tidal de- 
formations. In Section 5 we compute the light curves that 
arise from the numerical configurations (the numerical 
procedure used to compute the light curves is described 
in APPENDIX B). We explain how the reflection prop- 
erties of the materials affect the observed variability, and 
demonstrate the differences between our light curves and 
those resulting from the classical Roche approximations. 
Finally, in Section 6, we use our models to fit the ob- 
served properties of Kuiper Belt binary 2001 QG298, and 
to derive the physical properties and bulk density of this 
system. We summarize in Section 7. 

2. NUMERICAL METHOD 

We are interested in finding the self-consistent equi- 
librium configuration of tidally-locked homogeneous bi- 
naries, rotating in circular orbits. We assume that the 
bodies are strengthless, so that their equilibrium config- 
urations are determined by the conditions that the total 
(gravitational + rotational) potential is constant along 
the surface of each of the bodies (but may differ between 
the two bodies). 

The parameters that determine the equilibrium config- 
urations are the mass ratio q = M2/M1 < 1; the separa- 
tion between the components d; and the scaled angular 
velocity, uj 2 /Gp. Here we treat the separation and an- 
gular velocity as two independent parameters, and solve 
for the total mass of the system. For non-spherical con- 
figurations, the mass, separation, and angular velocity 
are not related through Kepler's law, and our choice of 
parameters allows us to conveniently explore the non- 
Keplerian nature of the solutions. Our numerical pro- 
cedure is described in detail in [X] We summarize the 
numerical scheme below. 

We parametrize the surface of each body by specifying 
the distance from its center to the surface along N preset 
angular directions, Ri=i...N. We associate each of the ./V 
surface-points with the small surface area surrounding 
it, and with a "mass-cone" stretching from the center 
to this surface area. The mass in the i'th cone is then 
Mi = AilpRf/3, where AO is the solid angle that the 
cone occupies. We distribute the N points evenly, so that 
each cone covers an equal solid angle as viewed from the 
center of the body. 

The potential at a point located on the surface is then 
the sum of the gravitational potential induced by the 
mass in all cones of both components, and the rotational 
potential about the common center of mass. For a model 
with Ni points sampling the surface of Mi , and N2 points 
sampling the surface of M2, our goal is to solve for the 
N\ + N2 values of Ri for which the total potential along 
the surface of each body is constant. The value of the 
potential may differ between the two bodies. 

We explicitly compute the gravitational and rotational 
potentials (seeEJ. Our algorithm solves for N\ + N 2 + 2 
variables, namely the values of the distances Ri defining 
the surfaces of the bodies, and the values of the constant 



potentials along their surfaces, ci on Mi and C2 on M2. 
The equations that we solve are the iVi + N2 conditions 
of constant potential, and the additional constrains pro- 
vided by the given mass ratio q, and separation d, 




E ieMl ,M 2 ^(%) + ^ot(fli) + c* , forj'EM,, 

M 1 /M 2 - q 

ZCoM.l - ICoM,2 + d 



(1) 

(see|A|, where a;coM,fc is the center of mass position of 
Mk, along the direction separating the two components. 

We use a Newton-Raphson scheme to solve for the vari- 
ables, x = (Rj,ci,C2) that satisfy equations (p}. We start 
with an initial guess for the configurations, and in each 
iteration correct the current guess by a small amount 
dx = — Fx J -1 , where J is the Jacobian derivatives ma- 
trix. We compute the Jacobian analytically, and iterate 
the Newton-Raphson scheme until the solution has con- 
verged so that F = to within some numerical threshold, 
and the correction dx is small compared to x. 

We define the x direction to be pointing from the more 
massive component to its companion, and the z direction 
to be parallel to the angular velocity vector. We assume 
that the equilibrium configurations have reflection sym- 
metries about the x—y and x—z planes. We therefore dis- 
tribute our N points on the surface of a quarter-sphere. 
In the discussion that follows, we use ~ 1600 points on 
a quarter-sphere for each of the components. 

3. EQUILIBRIUM FIGURES OF ROTATION 

3.1. Physical Review and Analytic Approximations 

In this section, we provide a brief introduction to the 
basic physical concepts underlying the derivation of equi- 
librium figures of rotations. We start by introducing the 
properties of single rotating Maclaurin spheroids. We 
then proceed to examine the properties of infinitesimal 
spheroids orbiting a massive companion. While the ex- 
act equilibrium configuration of a satellite may be far 
from spheroidal, this exercise captures the basic prop- 
erties of the equilibrium figures of binary systems. We 
then proceed to examine triaxial ellipsoids. We write 
the equations for the single Jacobi ellipsoids, and then 
study infinitesimal triaxial ellipsoids satellites in binary 
systems. 

In section f3. 21 we present our numerical results for the 
general equilibrium configuration of mutually deformed 
binaries. Our numerical configurations deviate from pure 
ellipsoidal forms. We discuss the properties of our so- 
lutions and compare them with the previous estimates 
based on the Roche ellipsoids in section |H 

3.1.1. Maclaurin Spheroids 

In his classical derivation of the rotational distortion 
of the Earth, Newton argued that since the Earth is in 
equilibrium, the weights of liquid filling two canals, one 
stretching from the center of the Earth to the equator, 
and the other from the center to the pole, must be equal. 
This is equivalent to demanding that the equator and 
the pole are a part of the same equipotential surface. 
Newton further pointed out that since the accelerations 
associated with both gravity and rotation are propor- 
tional to the distance from the center of the Earth, these 
weights (or, equivalently, the surface potentials measured 
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relative to the center) are given by 0.5 h a(h), where h 
is the height of the surface, and a(h) the acceleration on 
the surface. 

When Maclaurin later addressed the general problem 
of spheroidal equilibrium configurations, he followed the 
arguments outlined by Newton. For the potentials on 
the equator and at the pole we have, 

1 / 

1 9 equator 



equator 



-a X 
2 



u> 2 a) 



1 



•"pole — 2 C ^ 9pole 



(2a) 
(2b) 



where (/equator is the gravitational acceleration on the 
equator, g po ic is the gravitational acceleration at the 
pole, u) is the angular velocity, a the semi-major axis 
(a = b > c) of the spheroid, c = ayl — e 2 is the semi- 
minor axis, and e is the eccentricity. 

For a spheroid, the gravitational acceleration along the 
equator and at the pole are analytical functions of the ec- 
centricity e, and are given by (e.g. Chandrasekhar 1969), 



ffequator = 2nGpa 



sm e — e 



9 P olc = 4nGpa 



VT- 



sin e 



(3a) 



(3b) 



Requiring that the potentials on the equator and at the 
pole are equal then yields, 



2 _ 3cquator -flpolc Vl - E 2 

a 

2(3-2e 2 )sin- 1 e- £ (l - e 2 ) 



irGp 



(4) 



providing the relation between the eccentricity and an- 
gular velocity for Maclaurin spheroids. 

3.1.2. An Illustrative Example: An infinitesimal Spheroid 
in a Binary System 

Let us now consider the problem of a binary system, 
comprised of a massive primary of mass M orbited by 
an infinitesimal spheroidal (a = b > c) satellite. Here we 
assume that the system is tidally locked, so that the spin 
and orbital angular velocities are equal. For an infinites- 
imal satellite, the system's center of mass is effectively at 
the center of the massive primary, and we assume that 
the satellite is in a circular Keplerian orbit. 

The massive primary exerts a tidal force on the 
spheroidal satellite. On the point facing the primary, 
this results in an acceleration opposing the satellites' 
own gravitational acceleration, of mag nitude 2GMa/d 3 , 
where d is the separation between the two components. 
On the polar axis, the primary serves to compress the 
secondary, providing an acceleration along the polar axis 
of GMc/d 3 . 

Following the previous line of argument, the potentials 
on the surface of the secondary can be readily computed. 
Writing V a for the potential on the equatorial point fac- 
ing the primary, and V c for the polar potential, 



., 1 / 2 2GMa\ 
= 2~ a V 5oquator - w a d? — j 

1 / GMc 

Vc = r( 9 >° u + ^rd 



(5a) 
(5b) 



Demanding that these two potentials are equal, and 
assuming Keplerian rotation, 



GM 



we now obtain 



2 (/equator 
UJ = 



SpoicVl - e 2 



(4 - e 2 )a 



(6) 



(7) 



which relates the eccentricity and rotational velocity of 
infinitesimal spheroids in binary systems. We note that 
the resulting spheroids are not genuine equilibrium con- 
figurations, because even for the eccentricities dictated 
by equation 0, for which V a = V c , the value of the 
potential varies along the surface the spheroids. 
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Fig. 1. — The equatorial (solid) and polar (dashed) potentials 
versus eccentricity for an infinitesimal spheroid in a binary system. 
An equilibrium solution arises where the two curves cross (filled 
circles). The upper panel is for ui 2 /Gp = 0.377, the middle is for 
the "Roche limit" spheroid (see text) u) 2 /Gp = 0.4515, and the 
lower panel is for oj 2 /Gp = 0.503. In the shaded area (AV < 0), 
the net force is away from the center. 

Figure [T] shows the equatorial (V a ) and polar (V c ) po- 
tentials (relative to the satellite center, see equation [5]) , 
as functions of eccentricity for three representative val- 
ues of the angular velocity. The solid curves show the 
equatorial potential, and the dashed curves show the po- 
lar potential. Equilibrium configurations arise where the 
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two potentials obtain equal values and the curves cross. 
This is shown by the filled circles in Figure [T] 

For low angular velocities (panel a) there are two eccen- 
tricity solutions for a given angular velocity. The lower 
eccentricity solution is stable. A flattening of the equi- 
librium spheroid will result in restoring forces acting to 
raise the poles and compress the equator, and vice versa. 
The higher eccentricity solution is unstable. A flattening 
of the equilibrium spheroid will result in forces that act 
to flatten it even further, driving the eccentricity away 
from the equilibrium value. Figure Q] shows that the po- 
lar potential is always positive (relative to the satellite 
center), and the forces on the pole are therefore always 
directed toward the center. However, the equatorial po- 
tential on the point facing the primary becomes negative 
at large eccentricities, resulting in an outward force act- 
ing to unbind the material. 

Panel (b) shows that there is a limiting angular velocity 
for which only one solution exist. This is known as the 
"Roche limit" configuration. For even higher angular 
velocities (panel c), the potential at the pole is higher 
than the potential at the equator for any eccentricity, 
and no equilibrium solution can be found. 

Note that at the "Roche limit", the potential is still 
positive, and the material remains bound. The Roche 
limit is therefore crossed before the "Roche lobe"0 is 
filled. For larger values of angular velocity, the force on 
the pole is always larger than the force on the equator 
(see panel c), and the configuration is therefore driven 
to an ever increasing eccentricity. Eventually, the po- 
tential on the equator will become negative, resulting 
in mass shedding along the equator. Even though the 
Roche limit occurs before the Roche lobe is filled, once 
the Roche limit is crossed the configuration will evolve 
to shed mass. 

While our example assumed an unrealistic spheroidal 
configuration for the infinitesimal satellite, we find that 
it captures much of the physics at play in homogeneous 
binary systems. This includes the existence of two solu- 
tions, one of which is stable and the other not; the exis- 
tence of a "Roche limit" that is reached before the Roche 
lobe is filled; and the understanding that despite this 
fact, beyond the Roche limit mass shedding will even- 
tually occur. As we discuss below, the same principles 
apply for triaxial configurations. 



3.1.3. Triaxial Ellipsoids: The Jacobi and Roche Solutions 

Triaxial ellipsoids (a > b > c) are defined by two ec- 
centricities, ei = \J\ — (b/a) 2 and e 2 = \J\ — (c/a) 2 . 
To find the equilibrium configurations of single triaxial 
rotating ellipsoids, we demand that the surface potential 
along the three principal axes be equal. This provides 
two equations, which can be solved for the two eccentric- 
ities. 

The gravitational accelerations along the principle axes 
of a triaxial ellipsoid are (e.g. Chandrasekhar 1969), 



2irGpAiXi 



(8) 



4 The Roche lobe is the region of space around a body within 
which a test particle is gravitationally bound to that body. 



where 
Ai = a 1 a 2 a 3 



du 



o (a 2 + u)yj(a\ + u)(a\ + u)(a\ + u) 



(9) 

A comparison of the potential along the semi-major 
axis a, and the other equatorial axis 6, yields 



g a -up-a = g b \j\-e\-u?a(\-e\), (10) 
whereas for the minor (polar) axis c, 



9a 
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(11) 



These two equations can be solved numerically for the 
values of e\ and e 2 given an angular velocity uj. These 
solutions are the Jacobi ellipsoids. 

If an infinitesimal triaxial ellipsoid is rotating about a 
massive spherical primary, the potential must be modi- 
fied to include the tidal and gravitational contributions 
of the primary. In this case, 



g a — oj a 



WT^-o; 2 a(l- e ?) + ^Kl-e?) 



(12) 



and 



2 2GM a 
9a~ w a ^g— 



9c\/l-el + 



GMa 



(1-4). (13) 



Assuming Keplerian rotation (equation [S]), we may ex- 
press the above equations in terms of the eccentricities 
ei and e%, and the angular velocity only, 



9a - 9b 
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9a - g c yl - e 
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(14a) 
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These are the Roche ellipsoidal solutions for an extreme 
mass ratio (q <C 1). Figure [5] shows the Roche ellip- 
soid solutions for q <C 1. Panel (a) shows the angu- 
lar velocities given by equations 114a and b (dark and 
gray contours) as functions of the eccentricities e\ and 
e 2 . Equilibrium solutions exist when the two equations 
yield equal values of angular velocity. Examples of equi- 
librium configurations are shown by the filled symbols. 
The circles correspond to uj 2 /irGp = 0.06, and the trian- 
gle to uj 2 /nGp = 0.0901. 

Panel (b) shows the positions of the equilibrium so- 
lutions (solid curve). The dashed curve shows the line 
ei = e 2 (a prolate spheroid) for comparison. In panel (c) 
we display the angular velocity for which an equilibrium 
solution is found as a function of e±. As in the case of the 
spheroidal satellite, for low angular velocities, two equi- 
librium solution are found (see filled circles). Panel (c) 
shows that there exists a maximal value of angular veloc- 
ity for which only one solution can be found. This is the 
"Roche limit" ellipsoid, occurring at uj 2 /irGp = 0.0901 
(for q <C 1), shown here by the filled triangle. For values 
of angular velocity beyond the "Roche limit" , no solution 
can be found for any combination of eccentricities. 

3.2. Numerical Results 
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B 1 

Fig. 2. — Roche Ellipsoid solutions for q <C 1. Panel (a) shows the 
angular velocities given by equations 1141 (dark and gray contours) 
as functions of the eccentricities e\ and 62- A solution is found 
when dark and grey contours of equal values cross. For example, 
the filled circles show the two solution obtained for oj 2 /nGp = 0.06. 
The filled triangle is the single solution obtained for u> 2 /irGp = 
0.0901. The solid curve in panel (b) shows the eccentricities of 
the Roche ellipsoid solutions. The solutions shown in panel (a) 
are displayed here again. Panel (c) shows the angular velocity 
associated with these solutions. The single solution for u) 2 /nGp = 
0.0901 has the maximal possible angular velocity and is therefore 
the "Roche limit" ellipsoid. 

We have computed the equilibrium configurations of 
rotating binary systems for mass ratios, q, between 
10 -3 and 1. For each value of q, we found solutions 
starting at a very low angular velocity (large separation) , 




0.5 ■ < '''' 1 




D 



Fig. 3. — Equilibrium configurations of rotating binaries, (a) 
q = 1, distant components, (b) q = 1, close components, (c) 
q = 0.2, close components, (d) q = 0.01, close components. 

for which the rotational and tidal deformations are small, 
and following a sequence of increasing angular velocity, 
ending at the Roche limit where an equilibrium solution 
can no longer be found. 

In this section we present examples of our numerical re- 
sults for the equilibrium figures of rotation. These results 
were obtained using the numerical method described in 
section [2j which allows the bodies to take any form that 
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is symmetric about the x — y and x — z plane£0. As we 
demonstrate below, in some cases our numerically com- 
puted equilibrium figures depart from an ellipsoidal form. 

In figure we show several examples of equilibrium 
configurations. Panel (a) shows the equilibrium forms 
of equal-mass components, orbiting at a large separa- 
tion (small angular velocity). The equilibrium figures in 
this case are, as expected, close to spherical. Panel (b) 
shows the same system orbiting at a closer separation 
(note different scale). Asphericity becomes apparent for 
this rapidly rotating system. The equilibrium configu- 
rations resemble triaxial ellipsoids, but departures from 
pure ellipsoidal forms can be readily seen. 

Panel (c) shows an example of a smaller mass ratio, 
q = 0.2. In this case, the larger body is only mildly de- 
formed, due to the smaller gravitational perturbation of 
its low-mass companion. However the smaller body is 
visibly distorted by the gravitational forces of its high- 
mass companion. Departures from ellipticity are appar- 
ent for the lower-mass component. 

Finally, in panel (d) we show an extreme mass ratio, 
q = 0.01. For this system the center of mass is close to 
the center of the massive component. The lower mass 
satellite, is significantly deformed, yet remains nearly el- 
lipsoidal (see Section 4.1). The massive component is 
unaffected by the minute gravity of its low-mass com- 
panion, and is deformed merely by its own rotation, thus 
taking the form of the oblate Maclaurin spheroid of the 
appropriate angular velocity. 

4. DISCUSSION 

4.1. N on- ellipticity 

Consider the gravitational potential that a primary of 
mass M creates on the surface of its companion, which 
has a semi-major axis a, and is located at a distance d. 
On the axis connecting the two components, 

Tr , . GM T a fa\ 2 /a\ 3 , . 

^—i 1 -^) H- (15) 

The first term is constant, and does not produce any 
forces. The second terms gives rise to the circular orbital 
motion. The third term is the lowest order tidal force, 
which produces the symmetric deformations that were 
considered in the analysis of the Roche equilibrium ellip- 
soidal configurations (section |3.1.3[> . The fourth term is 
the first asymmetric correction to the tidal potential, and 
is a factor a/d smaller than the previous term. Asym- 
metric deformation are therefore always small at large 
separations. In addition, the primary is not spherical, 
as it is also deformed by rotation and by tidal interac- 
tions. This affects the gravitational potential, and de- 
viation from equation (|15l) may also contribute to non- 
ellipsoidal deformations. 

At a given separation, the semi-major axis of the 
lighter component is oc q 1 ^ . For small values of q, the 
asymmetric correction to the potential is therefore small, 
and the bodies remain ellipsoidal to a good approxima- 
tion. For larger values of q, the asymmetric correction 
become significant Figure 2] shows an example of the 

5 The bodies are separated along the x-axis, and the rotation is 
about the z-axis. 

6 For eccentric orbits, the average gravitational potential devi- 



departures from a symmetric ellipsoidal form, for q = 0.3, 
d = 2.636, and u 2 /Gp = 0.29732. The shaded area is the 
projection of our numerically computed equilibrium con- 
figurations on the x — z plane. The black contours show a 
projection of the best-fit ellipsoids. This example shows 
the asymmetric deformation of the lighter component. 
There is a significant "bump" on the side facing the pri- 
mary, while the distant side is somewhat depressed. The 
primary's asymmetric deformation is minor. 

To estimate the extent to which our equilibrium forms 
depart from pure ellipsoids, we fitted an ellipsoid to each 
numerically computed configuration, and then consid- 
ered the differences between the best-fit ellipsoid surface 
and the numerical solution. In figure [5] we display the 
departures from ellipsoidal configurations. The upper 
panels show the root mean square (RMS) of differences 
between the two forms, y/^2(Ri, c i — Ri,soi) 2 /N, where 
Ri : ci is the distance from the center to the best-fit ellip- 
soid surface, Ri, so \ is the distance from the center-of-mass 
to the surface of the numerical solution, and i = 1, .., N. 
The left hand panel is for the primary, and the right 
hand panel for the secondary. The RMS of the primary 
(secondary) and the separation were normalized to the 
primary's (secondary's) mean radius. In each panel, dif- 
ferent curves are shown for different values of q. The 
lowest curve is for the smallest mass ratio that we con- 
sider, 10~ 3 , followed by q — 0.01, 0.1, 0.2,..., 1.0 in 
increasing order. Some labels are shown near the curves 
for guidance. 

Figure [5] shows that indeed ellipsoids are a good fit at 
large separations, and become inaccurate at small sepa- 
rations. It also shows that for both of the components, 
the RMS at a given distance increases as a function of 
q. At small separations and comparable mass ratios 
(q > 0.2) the typical RMS is > 1%, and can be as large 
as ~ 3% for the most extreme cases. 

The departures from ellipsoidal forms are most promi- 
nent in the direction facing the companion (see Figure[4]). 
The RMS may therefore remain relatively small even in 
cases where locally, large deviations occur. In the lower 
panels of figure [5] we therefore show the largest devia- 
tion between the equilibrium configuration and its best 
ellipsoidal fit, for the primary (left panel) and secondary 
(right panel). At small separations and comparable mass 
ratios, deviation of order 10% (20%) occur for the pri- 
mary (secondary) component. 

4.2. Departures from Keplerian Rotation 

The tidal and rotational deformations of the equilib- 
rium configurations modify their gravitational potentials, 
causing them to deviate from those of spherical masses. 
These modified potentials affect the dynamics of the or- 
bit, leading to departures from Keplerian rotation. For 
a given total mass and angular velocity, the separation 
is generally larger than that predicted by Kepler's law. 

At large separations (small angular velocities), the 
bodies remain approximately spherical and the orbit is 
nearly Keplerian. At small separations (high u>), even for 

ates from a Keplerian potential, not only because of the deforma- 
tions of the body inducing the potential, but also because of the 
changing distance along the eccentric orbit. If this distance change 
is much smaller than the deformations, the circular approximation 
holds (e.g. for the case of equal mass, near-contact binaries this 
implies e < 0.01). Our models do not apply to higher eccentricities. 



8 



Gnat & Sari 



1 - 



N 




■1.5 -1 -0.5 0.5 1 1.5 2 25 3 3.5 

K 



Fig. 4. — Departures from ellipsoidal form, for q = 0.3, d = 2.636, ui 2 /Gp = 0.29732. The shaded area is the projection of the equilibrium 
configurations on the x — z plane. The black contours show a projection of the best-fit ellipsoids. 
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Fig. 5. — Departures from ellipsoidal configurations. The upper panels show the RMS of differences in radii between our numerical 
configurations and triaxial ellipsoids fitted to them, as a function of separation. The lower panels show the maximal local deviation 
between the numerical results and the best- fit ellipsoids. The left hand panels are for the massive component (primary), and the right hand 
panels are for the secondary. Different curves are for different mass ratios, q = 10" 3 , 10~ 2 , 0.1, 0.2, 0.3,... 1. Some labels arc indicated 
near the curves for guidance. 
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Fig. 6. — Departures (solid curves) from Keplerian rotation 
(dashed line) for mass ratios of 10 -3 , 10~ 2 , 0.1, 0.2, ... 1.0 (top to 
bottom). Some labels are indicated near the curves. 

extreme mass ratio (q <C 1), the rotational deformation 
of the primary leads to small departures from Keplerian 
rotation. For larger values of q the effect grows, as tri- 
axial and higher order deformations become significant. 

Figure [S] shows the extent to which our numerical solu- 
tions depart from Keplerian rotation, measured in terms 
of GM to t/u 2 d 3 . Keplerian rotation is given by the con- 
stant line GM to t/u 2 d 3 — 1 (dashed gray line). The nu- 
merical solutions for mass ratios between 10 -3 and 1 are 
shown by the solid curves. Departures from Keplerian 
rotation are a growing function of the angular velocities 
for any q. For q = 10~ 3 , the maximal departure near the 
Roche limit is of order 1%. For comparable masses, de- 
partures of order 10% occur at the Roche limit, reaching 
a maximal value of ~ 13% for equal mass components. 

4.3. Angular Momentum 

In figure [7] we show the angular velocity (ui 2 /AnGp) 

versus the angular momentum ( J/V47rGp 3 ^ 2 ^ 5 ^ 3 , where 
V is the total volume). Different curves are for different 
mass ratios, between 10 -3 (leftmost curve) and 1 (right- 
most curve). 

At a given mass ratio, the angular momentum first 
decreases as u) grows, but later begins to increase again. 
Each curve ends at a "Roche limit" appropriate for its 
mass ratio, where an equilibrium configuration can no 
longer be found. 

Focus first on the data for q = 10~ 3 , shown again in 
Figure |U For this extreme mass ratio, the system is 
rotating about the primary's center. Consider a primary 
of mass M and radius R, orbited by a satellite of mass 
m = qAI at a distance d. The angular momentum is the 
sum of two terms, an angular momentum due to primary 
spin, 

Jspin ^ ^MR 2 uj (16) 
and an angular momentum due to the secondary orbit, 
■/orbit - md 2 u>. (17) 
For Keplerian rotation, J ox bit oc w -1 / 3 . 
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Fig. 7. — Angular velocity versus angular momentum for our 
numerical solutions. Data shown for q between 0.001 and 1, as 
indicated by the labels. Lines are shown to guide the eye in con- 
necting points with the same q. 
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Fig. 8. — Angular velocity versus angular momentum. The solid 
and dashed lines show the analytic expressions for the spin angu- 
lar momentum of the primary, and orbit angular momentum of the 
secondary for 9«1. The thick gray line is the sum of both compo- 
nents. The points show the results of our numerical computations 
for q = 0.001. The inset shows the same data in linear scale. 

At low angular velocities, the separation between the 
two components is large, and the orbital angular momen- 
tum dominates. As the angular velocity increases, J or bit 
decreases, until the contribution of the primary's spin 
angular momentum begins to dominate, and the angular 
momentum is then proportional to uj. 

This is shown in figure [8] The dashed line shows the or- 
bital angular momentum (equation[T7|) that dominates at 
small angular velocities. The dark solid curve is the pri- 
mary spin angular momentum (equation ll€>[) that domi- 
nates at large w. The thick gray lines is the sum of both 
contributions. This simple analytic approximation nicely 
reproduces the trend followed by the numerical results. 

For less extreme mass ratios, the spin and orbit con- 
tributions due to both components must be taken into 
account. A simple analytical approximation may be de- 
rived under the simplifying assumptions of Keplerian ro- 
tation, and a known deformation. Here we assume that 
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the axes-ratios for both components of the binary sys- 
tem, and for both axes in each component are equal, 
such that bi/ax = c\ja\ — ^2/^2 = 22/02 = yl-r, 
where a, > hi > Ci are the principle axes of the i'th 
components. In this approximation, 

Jspin = ^Miof (1 + g 5/3 )(2 - e 2 )uj (18) 

and 

J orbit = &/ 3 M*' 3 q U -W (19) 
The total angular momentum is then given by 

J — Jspin + "/orbit- (20) 

At low angular velocities (large separations) the orbital 
angular momentum dominates, and J decreases with in- 
creasing to. However, as the separation decreases, the 
contribution of the spin angular momentum increases, 
and J grows again. 
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Fig. 9. — Angular velocity versus angular momentum for q = 1.0. 
The points show our numerical results. Th e different curves show 
the analytic approximation of equation 12(11 for various values of ec- 
centricity (e = 0,0.6,0.8,0.87,0.92), as indicated by the labels. As 
u increases tidal deformations grow, and the value of J continu- 
ously shifts to curves of higher eccentricity. 

For extreme mass ratios (q <C 1), the larger component 
dominates the spin angular momentum, and the smaller 
component dominates the orbit angular momentum. The 
asphericity is limited to small (spheroidal) rotational de- 
formations, and e remains close to 0. Equation I2TH thus 
agrees with its simplified version above (equations [16] and 
[T7|) . As the value of q grows, the tidal forces deform the 
bodies, and the configurations grow increasingly aspher- 
ical. 

In figure [5] we focus on the angular momentum for 
q = 1. The data points show our numerical solutions, 
and the different curves show the analytical approxima- 
tions of equation for q — 1 and for increasing values 
of eccentricity. The leftmost solid curve is for spheres 
(e = 0), and curves to the right are for higher eccentric- 
ities. At small angular velocities (large separations), the 
tidal and rotational deformations remain minor, and the 
spherical approximation agrees with the numerical re- 
sults. As oj increases, tidal forces deform the bodies into 



increasingly elongated shapes. Figure |H] indeed shows 
that as ui increases, the values of J grow further and fur- 
ther away from the spherical e = curve. In fact, as w 
increases the values of J continuously shifts to curves of 
higher and higher eccentricity. 

Returning now to figure [7J The angular velocity at 
which the angular momentum achieves its minimum is 
an increasing function of q. According to equation 1201 
this occurs at uj 2 cx q 3 / 2 (l + q)~ 2 (l +q 5 / 3 ) -3 / 2 . The an- 
gular momentum of equal mass binaries therefore "turns 
around" at a higher value of w than for binaries with 
extreme mass ratios. 

The approximation of equation thus explains the 
existence of a minimum to the angular momentum, the 
position of this minimum as a function of q, and the 
"drift" of the angular momentum to values larger than 
those expected for spheres. 

5. LIGHT CURVES 

We have carried out computations of the light curves 
of Kuiper Belt binaries, based on the equilibrium figures 
of rotation presented in Section IX2"1 The light curves de- 
pend on the relative positions of the observed binary, the 
observer, and the light source; on the inclination of the 
orbit relative to the plane of the sky; and on the reflect- 
ing properties of the materials composing the surface of 
the observed system. 

For binaries in the Kuiper Belt the observed geometry 
is simple. At a distance of > 40 AU, the observer (on 
Earth), light source (the Sun), and the observed Kuiper 
Belt binary are almost aligned, so that the "phase-angle" 
between the line-of-sight and the KBO-Sun direction al- 
ways remains small (< 2°). 

The reflecting properties of KBOs are not well con- 
strained. Here we consider two simple options. First, we 
consider uniform reflection, which produces an observed 
intensity proportional to the projected area on the plane 
of the sky. For the given geometry, this can result from 
simple "backscatter" reflectiorfl Second, we study the 
case of diffuse ( "Lambertian" ) reflection, for which light 
is reflected equally in all directions, and the observed in- 
tensity from a surface area depends only on the cosine of 
the angle between the Sun and the normal to the reflect- 
ing surface. 

The absolute intensity of light as seen from Earth also 
depends on the distances between the Sun, the KBO, 
and Earth, and on the KBO's size and albedo. Here we 
ignore the absolute intensities and instead focus on the 
normalized relative light-curves. We tile the surface of 
a given equilibrium configuration with surface triangles, 
and then compute the light intensity for a given orbital 
phase, inclination, and reflection law, taking into account 
possible obscurations of different tiles. The numerical 
procedure is described in detail in [B] We have verified 
that our results reproduce previous light-curve computa- 
tions given similar configurations (see Appendix IB. 2[) . 

In Figure [TUJ we show an example of a light curve 
computed for a mass ratio q = 0.6, an angular velocity 
uj 2 /Gp — 0.318, and an edge-on orbit. The solid curve 
shows the light curve for backscatter reflection, and the 

7 The Lommel-Seeliger law also reduces to being proportional to 
the geometrical cross-section at low phase angles (e.g. Lacerda & 
Jewitt 2007). 
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Fig. 10. — Light curves of an edge-on binary system with q = 0.6 
and uP- /Gp = 0.318. The solid curve is for backscatter reflection, 
and the dashed curve is for diffuse reflection. 



dashed curve is for diffuse (Lambertian) reflection. Both 
light curves have been normalized so that the intensity 
equals 1 at the maximum (note that for a given albedo 
the absolute intensity is larger for backscatter reflection). 

The light curve shape depends on the orbital parame- 
ters. In large separation (small u>) binaries, the eclipses 
only span a small fraction of the orbit, and the light 
curves generally appear to be slowly varying outside of 
the eclipses, as the projected area gradually changes. In 
close binaries (large to) , such "plateaus" do not exist be- 
tween the two eclipses. For the mass ratio considered in 
Figure ITUI uj 2 /Gp — 0.318 corresponds to a close binary, 
and the eclipses indeed appear to span most of the orbit. 

The mass ratio determines the depth of the eclipses. 
Even for the case of backscatter reflection, for which the 
intensity is proportional to the projected area, eclipses 
are deepest for equal mass components. In this case the 
tidal deformations are largest, and the projected area 
is a strong function of orbital phase. The deformations 
for extreme-mass ratio binaries are smaller (and more 
spheroidal), and the relative depth of the eclipses is thus 
smaller. For diffuse reflection the impact of the mass ra- 
tio is even larger. One of the most indicative differences 
between backscatter and diffuse reflections, is the rela- 
tive depth of the two minima. In backscatter reflection, 
which is sensitive only to the total projected area, the 
depths of the two minima are always equal, regardless of 
whether the small body is in front of the large body or 
vice versa. However for diffuse reflection, which is sen- 
sitive to the surface curvature, the two minima exhibit 
different depths. When the small body in viewed in front 
of the large body, a larger fraction of the surface is in- 
clined relative to the light source, and so the reflected 
intensity is smaller. An additional parameter that af- 
fects the depth of the eclipses is the orbital inclination. 
Figure [TU] shows the maximal variation that obtains for 
an edge-on orbit. Higher inclinations result in smaller 
variations. 

Next we examine the differences between light curves 
computed assuming the Roche approximation and light 
curves computed using our numerical solution. As an 
example, we consider the light curve of an equal mass 
binary, with an ang ular velocity to 2 /Gp = 0.316. In Fig- 
ure [IT] we show the ratio between the light intensity re- 



sulting from the Roche configuration appropriate to this 
angular velocity, and our non-ellipsoidal non-Keplerian 
configurations. Both light curves were computed using 
the numerical procedure described in Appendix B, as- 
suming edge-on inclination and diffuse reflection. They 
differ only by the equilibrium forms. The phase is de- 
fined so that the maxima of each individual light curve 
is obtained at and ir. 
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Fig. 11. — The ratio between the Roche-approximation and nu- 
merically computed-configurations light intensities as a function of 
orbital phase, for an equal mass binary with oj 2 /Gp = 0.316. In 
both cases we assumed an edge-on orbit. The solid curve is for 
diffuse reflection, and the dashed curve is for backscatter. The 
maxima of the individual light-curves occur at a phase of and n. 

Figurc[TTJshows that light curves differ by up to ~ 10%, 
and that the ratio between the two cases is a function of 
orbital phase. Intensity variations of ~ 10% can be easily 
detected with current facilities, and the numerical config- 
urations may thus produce superior fits to observed data, 
and provide better constrains on the physical parameters 
of the observed system. 

6. APPLICATION TO KUIPER-BELT BINARY 
2001 QG 298 

6.1. Observational Review 

2001 QG298 has been discovered as a part of the 
Hawaii Kuiper Belt Variability Project (Jewitt & Shep- 
pard 2002; Shcppard & Jewitt 2002, 2003). The ob- 
servations, carried out with the University of Hawaii 
2.2 m telescope and with the Keck 10 m telescope, are 
reported in Sheppard & Jewitt (2004; hereafter SJ04). 
SJ04 indicate that 2001 QG298 has a peak-to-peak pho- 
tometric range of 1.14 ± 0.01 mag in the ii-band, and 
a period of 6.8872 ± 0.0002 hours for a single- maximum 
period (which may arise due to albedo variations), or 
13. 7744± 0.0004 hr for a double-maximum period (which 
may arise due to rotation). SJ04 note that the double- 
peaked light curve provides a better fit to the observed 
data, for which the two minima differ by about 0.1 mag. 
Furthermore, Keck VBR colors of 2001 QG298 indicate 
no color variations along the period, within the photo- 
metric uncertainties of a few percent. 

The observed amplitude variation of 1.14 mag is ex- 
ceptionally large among large (> 25 km) solar-system 
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objects. SJ04 discuss three possible reasons for the large 
photometric range. First, they consider the possibility 
of albedo variations, and conclude that this scenario is 
unlikely. On asteroids, albedo brightness variations are 
usually smaller than 10— 20%. Larger albedo variations 
are likely to be associated with color variations (e.g. Ia- 
petus). Finally, the fact that a double-peaked period 
provides a better fit with two distinct minima favors a 
light curve produced by rotation rather than by albedo 
variations. 

Next, they consider the possibility of an elongated 
shape. SJ04 use the observed photometric range to in- 
fer an axis ratio b/a = 0.35 (assuming the intensity is 
proportional to the projected area). Given the abso- 
lute luminosity of 2001 QG298, they derive a semi-major 
axis between 170 and 270 km for "typical" KBO albe- 
dos range (0.04 — 0.1). Bodies of this size are unlikely 
to be held by material strength over long time-scales. 
Gravity-dominated bodies are rotationally deformed into 
Maclaurin spheroids or Jacobi ellipsoids, depending on 
their angular velocity. SJ04 note that the maximal pho- 
tometric range for stable Jacobi ellipsoids is 0.9 mag, 
lower than the observed variation of 2001 QG298- This 
maximum variation occurs for a Jacobi ellipsoid with 
b/a = 0.432 and c/a = 0.345 (Chandrasekhar 1969; 
Farinella et al. 1981) and more elongated objects are 
dynamically unstable. However, we note here that the 
maximum variation quoted in SJ04 applies to backscat- 
ter reflection, whereas for diffuse (Lambert) reflection 
these axis-ratios yield a variation of 1.5 mag. SJ04 also 
note that the rotational period of 2001 QG298 is too long 
to cause significant elongation for reasonable densities. 
They conclude that 2001 QG298 is unlikely to be a single 
rotating object. 

Finally, SJ04 consider the possibility of a close binary 
configuration. As discussed above, the components are 
distorted both by rotation and tidal forces. The pho- 
tometric range of 2001 QG298 is consistent with that 
of a comparable mass close binary Roche configuration 
(Leone et al. 1984). They conclude that given the large 
amplitude variation, long period, and difference between 
the two minima, a close binary is the most likely expla- 
nation for 2001 QG298- 

With the binary scenario, several attempts have been 
made to derive the physical properties of 2001 QG298 us- 
ing the Roche ellipsoidal approximations. SJ04 used the 
results presented in Leone et al. (1984) to estimate the 
density, and find p ~ 1 g cm~ 3 . Takahashi & Ip (2004) 
then constructed specific Roche binary light-curve sim- 
ulations to fit the observed light curves. Their best fit 
Roche solution implies that 2001 QG298 consists of two 
components with a mass ratio of 0.65. Their primary 
has axis-ratios b/a — 0.79 and c/a — 0.62, and their sec- 
ondary has b/a = 0.61 and c/a = 0.56. The separation is 
2.1 times th e prim ary's semi- major axis (see illustration 
in Appendix IB.2[) . This solution indicates a "mixed" re- 
flection pattern, with ~ 70% diffuse reflection and ~ 30% 
backscatter (uniform) reflection. For these parameters, 
they infer a bulk density of 0.63 ± 0.20 g cm' 3 . 

Later, Lacerda & Jewitt (2007) presented an additional 
model for 2001 QG298- In their solution the mass ratio 
is 0.84. Their primary has b/a = 0.62 and c/a = 0.65, 
and their secondary b/a — 0.45 and c/a = 0.41. The 



separation is 2.1 times the primary semi- major axis (see 
illustration in Appendix IB.2[) . Lacerda & Jewitt find 
that backscatter reflection best fits the observed data, 
and infer a bulk density of 0.590^Q 4y 3 g cm~ 3 . Lacerda 
& Jewitt also considered the possibility that 2001 QG298 
is in fact a single Jacobi ellipsoid with diffuse (Lambert) 
or backscatter reflection. They conclude that the Roche 
models fit the data significantly better than Jacobi ellip- 
soids. 

In the next section, we apply our numerical models 
to 2001 QG298, and demonstrate how our more accurate 
solutions can be used to constrain the physical properties 
of KBO binaries. 

6.2. Light-curve Fitting and Analysis 

We are interested in finding a physical configuration 
consistent with the observations of 2001 QG298- To this 
end, we compare the light curves associated with the 
equilibrium configurations computed in Section [3] with 
the observed light curve of 2001 QG298 (SJ04), after cor- 
recting for the light travel-time and phase angle effects. 

To find the best fit, we create a library of light curves 
for comparison with observations. However, because the 
light curve calculations are computationally expensive, 
we approached the fitting procedure in two steps. First, 
we created a library of "low-resolution configurations". 
Here we used equilibrium configurations computed with 
200 patches on a quarter-sphere, to construct light curves 
spanning the entire parameter space studied in Section|3J 
Given this library of light curves, we found the best fit 
(in the sense of minimum \ 2 ) model. 

For each model, which corresponds to a specific combi- 
nation of mass ratio, rotational velocity, and inclination, 
we fitted two parameters: a and /3. These are overall nor- 
malization factors for a "mixed" light curve composed of 
both backscatter and diffuse reflection, such that the to- 
tal intensity is 1 = aJ^f atter + ^J^f. In addition 
to these two normalization parameters, we allowed for a 
relative constant phase-shift between the observed and 
computed light curves. 




-0.5 ~'''-- ; < -D.5 



Fig. 12.— The best-fit configuration for 2001 QG 2 98- 

Given the best fit solution for the "low-resolution 
configurations" library, we computed a better-sampled, 
dense library of light curves based on "high-resolution 
configurations" , which have 1600 patches on a quarter 
sphere. For the high resolution configurations, we focus 
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on the parameter-space surrounding the "low-resolution" 
solution. We again searched for the best-fit solution 
among this new set of model light curves. 

The best fit solution is illustrated in Figure Q21 This 
model has a mass ratio, q = 0.93, and a rotational veloc- 
ity fl 2 /Gp = 0.333. The best fit is found for pure diffuse 
(Lambertian) reflection (so that a = 0), viewed at an in- 
clination of 3°. For the observed period of 2001 QG298, 
these parameters imply a density p = 0.72 g cm~ 3 . 
The observed magnitude then implies masses of 1.70 x 
10 18 (A/0.1)- 15 and 1.58 x 10 18 (A/0.1)- 15 kg for the 
two components, where A is the albedo. While our bod- 
ies are not ellipsoidal, an ellipsoidal fit yields semi-major 
axes of 102 x 78 x 71 and 102 x 75 x 69 km 3 . 

We display the observed light curve (symbols) along 
with the best fit model (solid curve) in Figure [TH For 
our assumed intensity errors of 3.5% on the computed 
light curves and 4% on the observed data, this models 
yields a x 2 of 123.4 for 107 degrees of freedonfl 




Fig. 13. — Observed light curve for QG298 (data points) and best 
fit solution (solid curve), obtained for q = 0.93, fl 2 /Gp = 0.333, 
with Lambertian reflection at an inclination of 3°. 



To estimate the errors on the best fit parameters, we 
inspect the values of \ 2 obtained for other models. With 
five fitted parameters and a minimal x 2 of 123.4, models 
with x 2 < 129.3 are within the la (68.27%) error region. 
Inspecting the model-library, we derive the following la 
error regions: For the mass ratio, 0.9 < q < 1.0; For 
the angular velocity, 0.332 < cu 2 /Gp < 0.334; And for 
the inclination, 2 < i < 4°. We further find that all 
acceptable models have "pure" diffuse reflection. Given 
the observed period (and error) , the error on the angular 
velocity can be used to compute the error on the bulk 
density. We find a 5.6% error on the density, p = 0.72 ± 
0.04 g cm -3 . 

7. SUMMARY 

In this paper we introduce a numerical method for 
computing the self-consistent equilibrium configurations 
of tidally-locked homogeneous binaries, rotating in circu- 
lar orbits. The equilibrium configurations depend on the 

8 For 112 observed data points and 5 fitted parameters, namely 
the mass ratio, angular velocity, inclination, values of a and 0, and 
the global phase shift. 



mass ratio, angular velocity, and separation between the 
two components. 

We explicitly compute the gravitational and rotational 
potentials on the surface of the two components, and 
use a Newton-Raphson-based scheme to converge to an 
equilibrium solution, for which the surfaces of each of 
the bodies is an equipotcntial surface. Our numerical 
procedure is described in detail in Section [2] and \K\ 

In Section [3] we discuss the properties of the equilib- 
rium figures of rotation. We begin with simple analytic 
approximations that we use to study the characteristics 
of the equilibrium solutions. We confirm that for low an- 
gular velocities (large separations) , two equilibrium con- 
figurations always exist, differing by their eccentricity. 
The low eccentricity solution is stable. A flattening of the 
low-eccentricity configuration results in restoring forces. 
The high-eccentricity solution is unstable. A flattening 
of the equilibrium configuration results in forces that act 
to flatten it even further, driving the eccentricities away 
from the equilibrium values. 

As the angular velocity increases, there exists a "lim- 
iting" angular velocity for which only one equilibrium 
solution exists. This is known as the "Roche limit" con- 
figuration. For even higher angular velocities, no equi- 
librium solution can be found. At the Roche limit, the 
material is still bound, and so the Roche limit is crossed 
before the Roche lobe is filled. However, once the Roche 
limit is crossed, the configuration is driven to an ever in- 
creasing eccentricity, and eventually mass-shedding will 
occur as material flows beyond the Roche-lobe. 

We find numerical solutions for the equilibrium config- 
urations for mass ratios, q, between 10~ 3 and 1. For each 
value of q we find solutions starting at a very low angular 
velocity (large separation), for which the rotational and 
tidal deformations are small, and following a sequence 
of increasing angular velocity, terminating at the Roche 
limit where an equilibrium configuration can no longer 
be found. 

Our numerical solutions indicate that the equilibrium 
configurations are not always ellipsoidal. Ellipsoidal fits 
are generally good approximations at large separations 
and for extreme mass ratios. Even at small separations 
and for comparable masses, the typical RMS deviation 
between a numerical solution and an ellipsoid fitted to it 
is 1 — 3%. However, departures from ellipsoidal forms are 
most pronounced in the direction facing the companion, 
and the maximal local deviations from ellipsoidal forms 
are of order 10% (20%) for the primary (secondary) com- 
ponent. 

The tidal and rotational deformations of the equilib- 
rium configurations modify their gravitational potentials. 
This, in turn, affects the dynamics of the orbit, leading to 
departures from Keplerian rotation. We measure depar- 
tures from Keplerian rotation in terms of GM to t/w 2 d 3 . 
At large separations the bodies remain approximately 
spherical, departures from Keplerian rotation remain 
small, and GM tot /w 2 (i 3 ~ 1 as it should. At small sepa- 
rations, departures from Keplerian rotation appear. For 
q = 10~ 3 the maximal deviation near the Roche limit is 
~ 1%. For comparable masses (q > 0.5), deviations of 
order 10% occur at the Roche limit, reaching a maximal 
value of ~ 13% for equal mass components. 

We inspect the angular momenta of the equilibrium 
configurations. We decompose the angular momenta into 
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"spin" and "orbit" components, and explain the trends 
in the angular momentum-angular velocity parameter- 
space. As the angular velocity increases, it shifts from 
being orbit-dominated to being spin-dominated, creating 
a minimum in the angular momentum. We discuss the 
increasing value of this minimum with increasing mass 
ratio, and demonstrate the "drift" of the angular mo- 
mentum to higher values as the tidal deformations grow. 

In section [5j we calculate the light curves that arise 
from our numerically computed equilibrium configura- 
tions, if placed in the Kuiper Belt. In addition to the 
configurations, the observed light curves also depend on 
the inclination of the orbit relative to the plane of the 
sky, and on the reflecting properties of the surface of the 
objects. 

We consider two possibilities for the reflecting proper- 
ties. First we consider an observed intensity that is pro- 
portional to the projected area on the plane of the sky, as 
is appropriate for backscatter reflection from Kuiper Belt 
objects. Second, we consider the possibility of diffuse 
(Lambert) reflection, in which the reflected light from a 
surface area is proportional to the cosine of the angle be- 
tween the sun and the normal to the surface. We note, 
that backscatter reflection is sensitive only to the total 
projected area, and therefore always yields two equal 
minima. However, diffuse reflection, which is sensitive 
to the surface curvature, generally produces two minima 
of different depths. When the small body is viewed "in 
front" of the large body, a larger fraction of the surface 
is inclined. 

We compare our light curves to those resulting from 
the classical Roche ellipsoidal approximations, and find 
phase-dependent intensity deviations of < 10% between 
the two cases. 

Finally, in Section |6] we apply our numerical models to 
Kuiper Belt binary 2001 QG298- This object exhibits an 



extremely large photometric range ofl.l4±0.01 mag (R- 
band), and a double-peaked period of 13.7744±0.0004hr. 
It is believed to be an example of a close-binary Kuiper 
Belt population (SJ04). 

Our numerical models indicate that 2001 QG298 is well 
fitted by a binary with a mass ratio q = 0.93^qq3, an 
angular velocity u 2 /Gp — 0.333 ± 0.001, a nearly edge- 
on inclination, i — 3 ± 1°, and pure diffuse reflection. 
For the observed period of 2001 QG298, these parameters 
imply a bulk density, p = 0.72 ± 0.04 g cm~ 3 

SJ04 estimate that 2001 QG298 is a representative of 
a large Kuiper Belt population of nearly-contact bina- 
ries, and that at least 10 — 20% of all large KBOs are, 
in fact, close-binaries. Upcoming LSST observations will 
identify > 20, 000 new KBOs (LSST Science Collabora- 
tions: Abell et al. 2009). If indeed a significant fraction 
of the large KBO population is in the form of contact 
binaries, the models and methods outlined in this paper 
may become essential for the interpretation of their light 
curves. 
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APPENDIX 

A. BASIC EQUATIONS AND NUMERICAL PROCEDURE 

A.l. Definitions 

We are interested in finding self consistent equilibrium configurations of a tidally locked binary systems, comprising 
of two homogeneous components, rotating in a circular orbits. 

The parameters governing the equilibrium configurations are the mass ratio q = M2/M1 < 1; the scaled angular 
velocity ui 2 /Gp; and the separation between the components center-of-masses, d. In Keplerian rotation, the angular 
velocity, separation, and total mass are related through equation f6|. However, here we allow departures from Kep- 
lerian rotation. When searching for equilibrium configurations, we treat the separation and angular velocity as two 
independent parameters, and solve for the total mass of the system. 

We first define a Cartesian coordinate system, whose origin is at the center of the more massive component, M\. We 
define the a; direction to point toward the lighter companion, M2, which is centered at x = d. The common center of 
mass is located at x = qd/(q + 1). We define the z direction to be parallel to the angular velocity. Figure IT41 illustrate 
the choice of coordinates. 

We describe the surface of each body relative to its own center of mass, using a spherical coordinate system. The 
surface of each component is sampled along N points with given angular directions, specified by the azimuthal angle 
<p and the polar angle 9 defined in the usual way. For example, in Figure [T4l 

x\ = ri sin(#i) cos((/>i) ( x 2 = d + r 2 sin(# 2 ) cos(^ 2 ) 

yi = n sm(0i)sm(<^i) I y 2 = r 2 sin(6> 2 ) sin(0 2 ) . (Al) 

z\ = r\ cos(6>i) { z 2 = r 2 cos(6> 2 ) 

The surface mapping is performed by associating each direction (4>i, 9i) with the distance, along this direction, between 
the center of mass and the surface, i?;. 

Each of the N points corresponds to a small surface area surrounding it. We distribute the N points evenly in <j) 
and in cos(#) to ensure that each point covers an equal solid angle Ml, as measured from the center of the body. We 
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Fig. 14. — A schematic representation of the coordinates used in the numerical solution. 



assume that the equilibrium configurations are symmetric about the x — y and x — z planes. This assumption allows 
us to only sample a quarter of a sphere when mapping the surface of each body. 

Each surface patch corresponds to a "mass-cone" stretching from the patch to the center of the body. For 
homogeneous bodies, the mass within the i'th "cone", stretching from r = to the surface where r = Ri, is 

Mi = /(f* pr 2 Andr = AflpRf/3. The total mass in the body is M = £\ Mi = AVLp/i x £\ Rf . 

The potential at a point located on the surface of a body is the sum of the gravitational potentials induced by the 
mass in all the cones in both components, and the rotational potential about the common center of mass. Our goal 
is to solve for the values of Ri, for which the total potential along the surface of each body is constant. While the 
surface of each body must be equipotential, the values of the potential may differ between the two bodies. 

A. 2. Evaluating the Potential 

To compute the gravitational potential at some surface point, we must first evaluate the distance between this point 
and a mass element dm inducing the potential. We can write the distance between two points that belong to Mi or 

M 2 , {(j>i,0i,ri) and (<pj,6j,rj) as 

Ar 2 = r\ +An + B (A2) 

with 

A = arj ± 2Zsin(6>;) cos(<^) 
B = l 2 + r 2 =F 2Wj sin(6>j) cos(^) 

where 

a = — 2{sin(#i) sm(0j)[cos((fo) cos(^) + sin(<^j) sin(<^,-)] + cos(#i) cos(0,-)} (A4) 

and 

. _ f d, i,j in different components (A5) 
— 1 0, i,j in same component ^ ' 

The plus (minus) sign in equation IA3al (|A3b[) should be used if j is in Mi and i is in M2, and the minus (plus) sign 
should be used if j is in M2 and i in M\. If both points are in the same mass, I = such that A = arj and B = 0. 
Given this distance, the gravitational potential due to the total mass in cone i with extent Ri, at a point (</>j, 9j, rj) 

is 

V i (r j )=AQ i {(f -M) ^/RjTARTTB + ^Vb 

+ - I ) ln (4 + ^ + VR't + + b) (A6) 

where A and B are the functions of rj specified by equation IA3[ and AQi is the solid angle covered by cone i. 



(A3a) 
(A3b) 
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For the rotational potential at a point ((f) j , 9j , rj ) we can write 

, 1 2 J (^) 2 +r J 2 sin(^) 2 -|^r J sin(^)co S (0 J ) , j in M x 

Krot - o w ^ / ( 2 • (A7) 

,iTTj + r j sin (^ ) 2 + ^ sin (^) cos (^) . 3 in M 2 

A. 3. The Newton-Raphson Scheme 

Consider a model with Ni points sampling the surface of Mi, and N2 points sampling the surface of M2. The 
conditions of constant-potentials provides N\ + N2 equation. For points Rj on the surface of M& these equations can 
be written as, 

ieMi i£M 2 

The variables in these equations are the iVi + N2 values Ri , and the values of the constant potentials on the surfaces 
of the two components, c\ on Mi, and C2 on M 2 . 

Two additional equations arc therefore required to solve the problem. One of these equations is the constrain 
provided by the given mass ratio, q, 



Afi 2 E ieMl R l 



A^i J2jeM 2 I'] 

r cones o 
idividual 

EieMi cos(0i) sin(^)i?f/4 E 3 eM 2 «»(&) Bin(^)J^/4 



-9 = 0, (A9) 



where Af^i and AO2 are the solid angles occupied by cones on Mi and M2, respectively. 
The last equation is for the distance between the individual center-of-masses of the two bodies, 



0. (A10) 



To summarize, the set of equation to be solved can be written as F = 0, where F is the N± + N2 + 2 dimensional 
vector, 

F jm (E t eM 1 ., 2 V l (Rj) + V Iot (R J ) + c 1 \ 

p(Nl+N2+2) _ _J F je ^ 2 II V i( R 3') + Kot(-Rj) + c 2 l=o. (All) 

f j Equation ( A9() for the mass ratio f 
(1) j Equation ( A10[) for the separation J 

The upper parentheses denote the size of each vector. 
In the numerical procedure, we start with an initial guess for the configurations, and iterate using a Newton-Raphson 

method (but see lA.5[) . correcting the values of our variables x — (Rj^Mn Rj£M 2 ) c ii c 2) by dx = —F x J -1 , where J is 
the Jacobian derivatives matrix. 

We compute the derivatives analytically. We iterate the Newton-Raphson scheme until the solution has converged, 
so that F = to within some numerical threshold, and the correction dx is small compared to x. 

A. 4. Symmetry 

As mentioned in Section IA.11 we assume that the equilibrium configurations are symmetric about the x — y and 
x — z planes. We therefore distribute N points on the surface of a 1/4-sphere, between = and n, and between 
cos(0) = 1 and 0. 

When computing the potential produced by the i'th cone with orientation (jU, = cos(0j), (f>i), the contributions of 
the three symmetric cones, at (/Zj, 2ir — 4>i), (— Mi; 4>i) an d (— Mi) 2tt — 4>i) must also be included. 
The vector F therefore becomes, 

F = J2i [ V Ri,^i,<t>ii R j) + VR ittH ,2n-(l>i(Rj) + V Ri,-Hi,<t>i{ R j) + VR u - fJ , i ,2Tr-<t> i (Rj)] + (A12) 
V I0 t (Rj ) + C = , 

and the Jacobian matrix J must be corrected accordingly to include the additional terms. 

A. 5. A Least Squares-Newton Raphson Solution 

We find that the numerical efficiency is much increased if instead of using equation IA10I for the distance between 
the center-of-masses, we use two equation, 

EigMi co8(&)sin(0i)i$/4 



EieM 2 R i/3 



(Al3a) 
0. (A13b) 
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We therefore solve an over-constrained problem, with Nl + N2 + 3 equations for iVi + + 2 variables. We substitute 
the Newton-Raphson matrix division dx = —F x J -1 , with a so-called "left-matrix division" dx — —F\3, which is 
the solution, in the least squares sense, to the over-constrained system of equations, J x dx = —F, We use the values 
of dx derived using this modified Newton-Raphson scheme to correct the current solution, and iterate to convergence. 

A. 6. Numerical Accuracy 

We consider the accuracy of our numerical method, as it applies to our best fit solution to the observed light curve 
of Kuiper belt binary 2001 QG 298 . This model has q = 0.93, n 2 /Gp = 0.333, and d = 2.8956. 

Given these parameters, we compare models constructed with various sampling resolutions, namely N = 
200,400,800,1200 and 1600 points on a quarter-sphere. We begin by estimating the "true" equilibrium configura- 
tions, by extrapolating the functions r^l/N) to (1/-/V) = for each direction i = 1 . . . N . For each N we then compute 

the RMS difference between the numerical solution and this "true" configuration, {Ri.so\ — Ri.truc) 2 /N, where 
Ri t sol are the distances from the center-of-mass to the surface in the numerical solution, and -R^true & r e the distances 
from the center-of-mass to the surface in the extrapolated, "true" solution. 

Since each solution is given at different azimuthal and polar angles, we first interpolate the center-to-surface distances 
of the high-resolution solutions onto the directions in which the lowest-resolution solution is known. We use these 
interpolated versions of the high-resolution solutions both to estimate the true configurations and to compute the RMS 
differences. 
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Fig. 15. — An example for the RMS differences between "true" equilibrium configuration, and numerical solutions at different resolutions 
(see text). The horizontal axis shows the number of points on a quarter sphere. The filled symbols are for the more massive primary, and 
the empty symbols are for the lower mass companion. 

Figure [T5] shows our results for the RMS difference versus the number of points. The filled symbols are for the 
more massive primary, and the empty symbols are for the lower mass companion. For example, the RMS difference 
between the primary solution with 1600 points and the "true" solution is 0.0048. We note, that the interpolation 
process introduces some errors that may increases this RMS. These values can be compared with the RMS difference 
between the numerical solution and the best fit ellipsoid (see Section 4.1 and Figure [5]) which are 0.023 and 0.027 for 
the primary and companion, respectively. 

B. COMPUTING LIGHT CURVES 
B.l. Numerical Procedure 

To compute the observed luminosity, we tile the surface of a given binary-configuration with triangles, typically two 
triangles for each of the N points used when computing the equilibrium configuring. For the computations presented 
in sections [3] and 01 which contain ~ 1600-point on a quarter-sphere, this corresponds to ~ 13,000 triangular tiles on 
each of the components. 

We compute the x, y, z coordinates of the vertices of the triangular tiles. Given the orbital phase and inclination, 
we then rotate the coordinate system such that the positive z-direction points toward the observer. Tiles with higher 
values of z are closer to the observer, and may block tiles with lower values of z from sight. 

We now sort the tiles of both bodies according to their mean z values. For each tile, we compute the area projected 
on the x — y plane, A*~ v , and the angle between the normal to the tile surface and the Sun, 

We define a "visibility" function for each tile, V$, describing the fraction of the tile's projected surface that is visible 
to the observer. For each tile, we estimate the visibility function by searching for other tiles which may block it 
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from sight. We start with the tile closest to the observer (highest z value), and search for possible overlaps between 
its x — y position and the x — y positions of all the lower- z tiles. When an overlap exists, we reduce the visibility 
function of the lower tile, by an amount corresponding to the fraction of its area that overlaps with the higher tile, 
Vi =Vi — ^4°™ llap /A*~ v . We then repeat this procedure for every other tile, correcting the visibility of lower tiles. 

This provides an overestimate of the reduction in visibility, since several high tiles may cover the same part of a 
lower tile, resulting in a "double-reduction" of the same covered area. This error remains small due to the simple 
shapes considered herfl 

The total observed light is given by ^ i A^~ y x Vi for the case of backscatter reflection, and by J^i ^i~ V x V{ x cos(i?j) 
for the case of diffuse (Lambertian) reflection. To compute the light curve due to a given configuration viewed at a 
given inclination, this process must be repeated for a set of orbital phases between and 2n. 

B.2. Sanity Checks 

To verify the accuracy of our light-curve computations, we compared our models against published results for similar 
configurations. In particular we focused on previous suggested models for QG298- Two previous published models 
exist, both using the Roche ellipsoidal solutions. The first is by Takahashi & Ip (2004), and the second by Lacerda & 
Jewitt (2007). 

According to Takahashi & Ip (2004), QG298 consists of two components with a mass ratio, q — 0.65. The primary 
has ei = 0.61 and ei = 0.69, and a secondary has ei = 0.79 and e2 = 0.83. The separation is 2.1 times the primary's 
semi-major axis. The system is illustrated in Figure [TBI 




Fig. 16.— The configuration found by Takahashi & Ip (2004) for QG 2 98- 

Takahashi & Ip find that their best fit indicates a ~ 70% diffuse reflection plus ~ 30% backscatter (uniform) 
reflection. For these parameters, they infer a bulk density of 0.63 ± 0.20 g cm~ 3 . 




0.4 0.6 
phase 

Fig. 17. — The light curve that we compute for the parameters inferred by Takahashi & Ip (2004). Compare to their Figure 2. 



Figure [T71 shows the light curve that we compute given the parameters inferred by Takahashi & Ip, along with the 
observed data points (left). A comparison of this light curve with that presented in Figure 2 of Takahashi & Ip (right), 
shows that our results are nicely consistent. 

According to Lacerda & Jewitt (2007), QG298 is a binary with a mass ratio q — 0.84. Their primary has e\ — 0.69 
and ei — 0.76, and their secondary e\ = 0.89 and ei — 0.91. The separation is 2.1 times the primary semi-major axis. 

9 We verified that we obtain the correct projected area for the com onent 

simple case of backscatter reflection from two spherical or elliptical 



Equilibrium Figures of Binary KBOs 19 

This system is illustrated in Figure [T5] Note that Lacerda & Jewitt accepted this configuration even though the two 
components overlap. They argue that they chose to accept it, because of the inaccuracies imposed by the ellipsoidal 
and Keplerian nature of the Roche approximation, and due to our poor understanding of the formation mechanisms 
of such binaries. 



Fig. 18. — The configuration found by Lacerda k. Jewitt (2007) for QG298- 



Lacerda & Jewitt find that backscatter reflection best fits the observed light curve. These parameters indicate a 
bulk density of O^O^ " 3 g cm 3 . In Figure [19] we display the light curve that we compute based on the parameters 
inferred by Lacerda & Jewitt (2007). Our light curves are again nicely consistent with those computed by Lacerda & 
Jewitt and shown in their Figure 9 (lower-left panel). 
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Fig. 19. — The light curve that we compute for the parameters inferred by Lacerda & Jewitt (2007). Compare to their Figure 9 (lower-left 
case). Here the solid curves are for backscatter reflection, and the dashed curves are for diffuse reflection. 
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